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ABSTRACT 

The thermal Sunyaev Zeldovich (SZ) effect directly probes the thermal energy of the 
universe. Its precision modeling and future high accuracy measurements will provide a 
powerful way to constrain the thermal history of the universe. In this paper, we focus 
on the precision modeling of the gas density weighted temperature Tg and the mean 
SZ Compton y parameter. We run high resolution adiabatic hydro simulations adopt- 
ing the WMAP cosmology to study the intergalactic medium (IGM) temperature and 
density distribution. To quantify possible simulation limitations, we run n = —1,-2 
self similar simulations. Our analytical model on Tg is based on energy conservation 
and matter clustering and has no free parameter. Combining both simulations and 
analytical models thus provides the precision modeling of Tg and y. We find that the 
simulated temperature probability distribution function and Tg shows good conver- 
gence. For the WMAP cosmology, our highest resolution simulation (1024'^ cells, 100 
Mpc/h box size) reliably simulates Tg with better than 10% accuracy for z > 0.5. 
Toward z = 0, the simulation mass resolution effect becomes stronger and causes 
the simulated Tg to be slightly underestimated (At ^ = 0, ~ 20% underestimated). 
Since y is mainly contributed by IGM at z > 0.5, such simulation effect on y is no 
larger than ~ 10%. Furthermore, our analytical model is capable of correcting this 
artifact. It passe_s all tests of self similar simulations and WMAP simulations and is 
able to predict Tg and y to several percent accuracy. For low matter density ACDM 
cosmology, the present Tg is 0.32(c^8/0■84)3■05-0•l5^^- (f7„/0.268)i-28-0-2'^B ^eV, which 
accounts for 10~^ of the critical cosmological density and 0.024% of the CMB energy. 
The mean y parameter is 2.6x 10"S((T8/0.84)-* i-2f^™(a„/0.268)i-28-0-2<^s. The current 
upper limit of y < 1.5 x 10~^ measured by FIRAS has already ruled out combinations 
of high (Tg > 1.1 and high flm > 0.5. 

Key words: Cosmic microwave background-theory-simulation: large scale structure, 
intergalactic medium, intracluster gas, cosmology, thermal history 



1 INTRODUCTION 

Ionized electrons with thermal motion can scatter CMB 
photons to generate secondary CMB temperature fluctua- 
tions known as the thermal Sunyaev Zeldovich (SZ) effect. 
Since all free electrons participate in the inverse Comp- 
ton scattering and contribute to the SZ effect, the SZ ef- 
fect is an unbiased probe of the thermal energy of the uni- 
verse at z < 6, for which the universe is highly ionized. 
The thermal SZ effect is sensitive to various physical pro- 
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cesses like adiabatic gravita tional heating, feedbac k, pre- 
heating and radiative cooling l|da S ilva et al.^2001(:|Lin et alJ 
120021: IWhite. Hernauist fc SpingeJl2002l : IZhang fc Wj2003h 
which affect the thermal energy of the baryons. In ad- 
dition, Compton cooling of first st ar supernova remnants 
(fOh. Coor ay fc KamionkowskillioO^) . cluster magnetic field 
tZhana 20041) . etc. could further alter the thermal energy of 
the universe to the level of > 10%. Therefore, the precision 
measurement and interpretation are of great importance to 
understand the thermal history of the universe. 

Current CMB expe riments such as CBIl B ond et all 
120021: iMason et al]l2003l). BIMA iPawson et alJl2002^ and 
ACBAR JKuo et al. 112003) marginally detected the SZ ef- 
fect. Several upcoming CMB experiments such as ACT, 
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APEX, Planck, SPT and SZA are likely able to measure the 
SZ effect with ~ 1% accuracy in the next several years. In 
order to utilize the power of such accurate experiments, the 
modeling of the SZ effect is required to meet this ~ 1% accu- 
racy. The first natural step for such modeling is to robustly 
understand the evolution of the baryon thermal energy in 
an adiabatically evolving universe. It is not only required 
to extract more complicated physics by comparing with ob- 
servations but also provides clues for the modeling of these 
complicated physics. 



Much effort has been devoted towa rd this 
goal, both anal yt ically jCole fc Kaiser| 



198S; 



199C; 



Makino fc Sutd Il99^ lAtrio-Barandela fc Mucket| 
Komatsu fc Kitavama '19991; 'Coorav, Hu fc Tegmark' 
20001: iMolnar fc B ir kinshaw. .2000 : Maiumdar JOOl; 
Zhang fc Penl T200ll: iKomatsu fc Sel iak 200 ^) an d sim- 
ulationally Jda Silva et al.l l2000l:~[Refregier et alJ l2000l: 
[ Seliak. Burwell fc PenI I2001I: ISpringel. White fc Hernauisd 
I2OOII: IZhang. Pen fc Wand l2002ll . Both methods have 
limitations, which has not been quantified and corrected 
to meet the precision of future observations. Analytical 
models of the SZ effect are often ad hoc procedures. In 
the halo model, the cluster gas pressure distribution is 
a free function of cluster mass and redshift. Though it 
can be calculated by various assumptions such as hydro- 
static equilibrium, its uncer tainty is hard to qu antify. In 
the continuum field model iZhang fc Penlliooj) . the gas 
temperature is determined by the gravitational potential, 
whose zero point is determined somewhat arbitrarily. So, 
analytical models must b e tested and calibra t ed ag ainst 
simulations. For instance, iRefregier fc TevssieI^ l|2002^ has 
tested the halo model against simulations and found good 
agreement. However, the conclusions drawn from these 
comparisons should be viewed with some caution. Numeri- 
cal simulations are known to have artifacts, stemming from 
limited resolution and finite volume. The impact of some 
of the artifacts have been investigated for the thermal SZ 
effect l^fregicr fc T cvssicr 2002) and the kinetic SZ effect 
JZhang. Pen fc Trajl20o3) . If not corrected, such artifacts 
would lead to biased calibrated analytical models. 



The SZ mean temperature decrement, or equivalently, 
the mean SZ Compton y parameter, which corresponds to 
the density weighted gas mean temperature Tg, are the low- 
est order SZ statistics. They are also the easiest to simulate 
and model. So, the precision prediction of Tg and y stands as 
the first natural step toward the precision modeling of the 
SZ effect. Their precision modeling also provides clues for 
the next low order SZ statistics such as the SZ power spec- 
trum and the corresponding gas pressure power spectrum. 
In this paper, we present a detailed study of the IGM den- 
sity and temperature distribution from a series of ACDM 
and self similar simulations. We further test the continuum 
model prediction of Tg and y parameter against simulations. 
Our goal is to quantify and correct numerical limitations and 
build calibrated analytical model aimed at 1% accuracy. We 
will follow a similar procedure as in this paper to discuss the 
precision modeling of the SZ power spectrum in a companion 
paper (Zhang, Pen fc Trac, 2004, in preparation). 



2 THE THERMAL SUNYAEV-ZELDOVICH 
EFFECT 

Free electrons scatter off CMB photons by their thermal 
motions and introduce secondary CMB temperature fluctu- 
ations; 



e = -2yST{u) = -2y 
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where x = /li^/fcsTcMB ■ Thi s is known as the therma l 
Sunyaev-Zeldovich (SZ) effect dZeldovich fc Sunvaeviri969t) . 
The Compton y parameter is given by the integral of elec- 
tron thermal pressure energy along the line of sight 
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where x is the comoving distance and a is the scale fac- 
tor. The lowest order statistics of the SZ effect is the mean 
Compton y parameter: 
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where Tg = {{1 + Sg)Tg) is the gas density weighted mean 
temperature and x = x/{c/ Hq) is the dimensionless comov- 
ing distance while Hq is the present Hubble constant. 

The thermal energy of the universe only accounts for a 
tiny fraction of the total energy of the universe. 



f^TE = 1.08 X 10" 



(4) 



O.SkeV 0.02 

As a comparison, f2cMB = 2.48 x 10~^h~^ and the energy 
in other wavelength bands of light 



f^EBL = 2.48 X 10"®/i"^ 



100 nw m~2 sr^i 



(5) 



3 HYDRO SIMULATIONS 

We ran cosmological hydrodynamical simul ations using 
a new Eulerian cosmological hydro code (|Trac fc PenI 
l2003allbl) . This Eulerian code (hereafter TP) is based on the 
finite-volume, fiux-conservative total variation diminishing 
(TVD) scheme that provides high-order accuracy and high- 
resolution capturing of shocks. The hydrodynamics of the 
gas is simulated by solving the Euler system of conservation 
equations for mass, momentum, and energy on a fixed Carte- 
sian grid. The gravitational evolution of the dark matter is 
simulated using a cloud-in-c ell particle-mesh (PM) scheme 
jHocknev fc Eastwoodlll988^ . 

The robustness of the TP code has been tested by com- 
paring the evolution of the dark matter and gas density 
p ower spectra from th e simulations with the fitting formula 
of ISmith et al. I ll2003f) . We also performed a code compari- 
son b y running the same initial conditions using the MMH 
code iPenlll99a) . which combines the shock capturing abil- 
ities of Eulerian schemes with the high dynamic range in 
density achieved by Lagrangian schemes. Power spectra are 
computed using FFTs. We find good agreement at all rele- 
vant scales and redshifts for both comparisons. 

Eulerian schemes are ideal for simulating the evolution 
of the IGM to model the thermal and kinetic SZ effects and 
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Figure 1. The gas temperature distribution in the 1024^, 100 
Mpc/h WMAP simulation. The top left panel shows the cumula- 
tive mass fraction p(> T) with temperature higher than T. The 
bottom bottom panel shows the conditional p(> T) in different 
overdensity regions at z = 0. The right panels plot the corre- 
sponding dP/d InT. At z = 0.0, only 5% gas is hotter than 1 keV. 
Gas temperature is tightly correlated with gas density. Nearly all 
virialized gas is hotter than 1 keV while effectively no gas with 
(5 < 10 is hotter than 0.1 keV. 

the Lyman alpha forest, because of their high speed, supe- 
rior mass resolution, shock-capturing abilities. Furthermore, 
Eulerian algorithms are computationally very fast and mem- 
ory friendly, allowing one to optimally use available compu- 
tational resources. 

We ran a 1024^ cells, 100 Mpc/h box size simulation 
with the best fit WMAP-alone cosmology Qm = 0.268, 

= 0.752 Oj, = 0.044, h = 0.71, and erg = 0.84 
JSoergel et alJl2003D . The ratio of dark matter particles to 
fluid elements is 1:8. We achieve a spacial resolution of 
Ax ~ 100 kpc//i and a dark matter particle mass reso- 
lution of Am ~ 5.6 X 10* Mq. The initial conditions are 
generated by sampling from an initial power spectrum com- 
puted using CMBFAST (Scliak & Zaldarriaga 1996). Tfiis 
simulation is started at a redshift of 2: = 100 and evolved 
down to z — 0, with data outputs at z = 3, 1, 0.5, 0.2 and 
1. This simulation takes approximately 700 time-steps to 
evolve from z — 100 down to 2: = 0. On a GS320 Compaq 
Alpha server with 32 cpus and total theoretical peak speed 
of 32 Gflops, the run requires 40 G memeory and takes 
approximately two days. Simulations are limited by both the 
box size, which causes the absence of large scale density fluc- 
tuation at scales larger than half box size, and resolution, 
which results in the failure to resolve small scale structures. 
Self similar simulations are ideal to test and quantify such 
simulation limitations since different redshifts directly cor- 
responds to different resolution and box size. We ran one 
n — —2 and one n = — 1 (fim — 1) self similar simulation 



Figure 2. The contour of gas (log)temperature and (log)density 
in the 1024^ WMAP simulation. Gas temperature strongly cor- 
relates with gas densit y with a sca l ing re l ation T oc p , as fo und 
in previous works (e.g. iKang et all il994h : fDave et al.l J200 Jl V 



with 512^ cells and the same amount of dark matter par- 
ticles. In these simulations, Qb is set to be 0.044/0.268 to 
mimic the Qt /^m ratio of the WMAP cosmology. The initial 
fluctuation is normalized such that, when linearly extrapo- 
lated to z — 0, the correlation length is half the simulation 
box size. We also run a 512^ cells, 100 Mpc/h WMAP sim- 
ulation with the identical initial condition as the 1024'^ for 
a direct comparison. The moving frame of the TP code is 
not turned on in these runs. We will check its effect in the 
future. 

We show the temperature distribution function, namely 
the mass fraction of gas hotter than T, p{> T) of the 1024"^ 
WMAP simulation in Fig. [H At 2 = 0.0, only 5% gas are 
hotter than 1 keV and ~ 40% gas are hotter than 0.1 keV. 
These fraction drops to 3% and 24% at z = 1.0. Gas tem- 
perature shows a strong positive correlation with its den- 
sity. At z = 0.0, nearly all gas with S > 1000 are hotter 
than 1 keV. For 100 < 5 < 1000, almost all gas are hot- 
ter than 0.1 keV. Nearly no gas hotter than 0.1 keV lies 
in (5 < 10 region. It is interesting to see how much gas lies 
in virialized halos. Virialized gas should have an overden- 
sity larger than that of the gas overdensity at the virial ra- 
dius. For an isothermal density profile p oc r~^, this states 
S > S(r^ir) = l/SAJQm ~ 100. ~ 100 for WMAP 
jEke et aLllljg^) is the mean matter density in a virialized 
halo with respect to the critical density. This factor 1 /3 does 
not change much for other profiles such as NFW and thus we 
omit its variation. Then the bottom panel of Fig. Q implies 
that only ~ 28% gas resides in virialized halos. At z — 1.0, 
this fraction drops to ~ 19%. 

In our simulation, a large fraction of gas is colder than 
~ lO^K. At 2; = 0, this fraction is ~ 10% and at 2 = 3, this 
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Figure 3. The comparison of the conditional temperature distri- 
bution function p(> T,Si < S < S2) between WMAP 1024^ (thick 
lines) and 512'^ (thin lines) simulations. For clarity, we only show 
S > 1000 (solid lines), 1000 > <5 > 100 (dot lines), 100 > (5 > 10 
(short dash lines), 10 > 5 > 1 (long dash lines). Though the over- 
all p(> T) agrees very well for two simulations (Fig.|5J, the condi- 
tional p{> T,Si < 5 < 62) does not converge yet and reflects the 
effect of simulation resolution. 5 > 1000 corresponds to the inner 
regions of clusters and groups. For these virialized objects, their 
temperature scales with respect to their mass M as M2/3. For 
clusters with virial radius around 1 Mpc/h, T ~ 1 keV. Higher 
resolution simulations are able to resolve more low mass halos, 
which are smaller and have lower T and thus result in a higher 
p(> T,5 > 1000) at T < 1 keV. The limited simulation resolution 
introduces artifacts into the T-p correlation. Though its effect to 
the gas density weighted temperature is minor (Fig.|S), it would 
affect the simulated pressure power spectrum significantly. 



fraction reaches ~ 50%. In reality, part of these gas may 
condense into stars or interstellar medium. Part of them 
may be photoionization-heated to above lO^K. Our simula- 
tion does not include any photoionization, radiative cooling, 
etc., so the prediction about these gas i s not reliable and i s 
hard to c ompare with other works (e.g. iKang et al.l Jl994l : 
IPave et~a l. (2001)). But the contribution of such gas to the 
SZ effect is negligible due to their low temperature, so we 
omit the complexity caused by such gas. 

The strong correlation of T and p observed in p(> T, > 
p) is clearly shown in the T-p contour of our 1024^ simulation 
(Fig. |5J. T oc p ho lds in a large T-p regions, as found in 
previous works (e.g. lKane et al.l l(l993 l: lDave et al.l l)2n0l|) ). 

Simulations are both mass and spacial resolution 
limited. So the above results may be strongly reso- 
lution d ependent. Both the Press-Schechter formalism 
Pr^^^ ^ Schcchtei||l974l and the Jenkins fitting formula 



Jenkin^t al.l 20^1) imply the existence of numerous small 



halos. The failure to resolve these small halos will result in 
an underestimation of the fraction of virialized gas and bi- 
ased T-p relation. To estimate the resolution effect, we com- 



Figure 4. The contour of gas (log)temperature and (log)dcnsity 
in the n = — 1 self similar simulation. The simulations show good 
convergence in most T-p regions. Higher resolution is required to 
simulate both low T, low p regions and high T, high p regions. 
But low T, low p regions have little contribution to mean gas 
density weighted temperature. 



pare between WMAP simulations and between self similar 
simulations. 

In low density regions, the density field is well resolved. 
But such regions generally have high Mach number and thus 
shocks are relatively poorly resolved. Since gas thermal en- 
ergy is generated by shock heating, in low density regions, 
temperature field is likely poorly resolved. In high density 
regions, it is the opposite case. Fig.|H]shows that both fields 
are indeed significantly affected by resolutions. The fraction 
of S > 1000 gas increases from ~ 1% in the 512^ simulation 
to ~ 5% at 1024^. The 1024^ simulation is able to resolve 
smaller halos, which have lower temperature and results in 
an increase in dp{> T,S > 1000)/dlnT at T < 3 keV. The 
fraction of virialized gas {5 > 100) increases from ~ 18% to 
~ 28%. For 5 < 10, the conditional p(> T) agrees well at 
the high T tail where Mach numbers are low and diverges 
at the low T tail where Mach numbers are high. So, for 

5 < 10, shock capture is the dominant resolution factor. In 
summary, the IGM density and temperature distribution in 
simulations is mainly limited by density resolution in high 
density regions and shock capturing ability in low density re- 
gions. We thus expect that higher resolution simulations will 
result in larger dispersion in the gas p-T phase space distri- 
bution, and thus larger pressure dispersion. This conclusion 
is further confirmed in the T-p contours of the probability 
distribution for our self similar simulation results (Fig. 
Such resolution effect has only minor effect on T, but its 
effect on the pressure power spectrum may be important. 
We will study this issue in a companion paper (Zhang, Pen 

6 Trac, 2004, in preparation). 

Our simulations suggest that (1) only a small 
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Figure 5. The comparison of the overall temperature distribu- 
tion function p(> T) between the WMAP 1024^ and the 512^ 
simulations. Two simulations agree well. p(> T) is directly re- 
lated to the mean gas mass weighted temperature Tg through the 
relation Tg = J Tdp. The good convergence of p(> T) implies 
that the simulation effects to Tg is minor, though we will still 
quantify them in Jl] 

fraction of gas are virialized and (2) a strong corre- 
lation exists between gas temperature and density. 
Both apparently contradict with the halo model. 
The halo model assumes that all gas resides in virial- 
ized halos. This assumption seems to contradict with 
(1), as noticed by Refreeier & Tevssier (200:^. The 
simplest version of the halo model assumes gas to be 
isothermal with temperature equal to its virial tem- 
perature. Since halos roughly have the same density, 
the predicted T-p correlation should be very weak. 
The weak dependence of halo concentration number 
on halo mass can not alter this straight prediction. 
More complicated int racluster gas models su ch as 
the universal gas profile dKomatsu fc: Seliakll2'00lf) pre- 
dict a weak variation of gas temperature from the 
core to the virial radius. But such variation is still 
too weak to explain the simulated T oa p relation. 

But simulations are resolution limited and some 
artifacts can be clearly seen in Fig. 3. It is likely 
that the part of the contradictions discussed above 
are caused by the limited simulation resolution. The 
limited simulation resolution effectively smooths the 
density field and may cause an underestimation of 
the fraction of virialized gas. The spacial resolution 
(~ 0.1 Mpc/h for our largest simulation) instead of 
the mass resolution is the main limiting factor to 
resolve halos. If our simulation only resolves halos 
with virial radius r^ir ^ rmin ~ 0.5 Mpc/h, following 
the Press-Schechter formalism (Press & Schechteg 
Il974h . such halos contain fh = l- Erf[-i^/V2] ~ 30% 



of the total mass at 2 = 0. Here, 1/ = Sc/o-{mmin), 
where Sc ~ 1.686 is the density threshold, a{m) is the 
density fluctuation in a sphere with mean mass m 
and m,nin is the mass of halos with rvir = ^min. This 
mass fraction is consistent with our simulation re- 
sults. But this predicted mass fraction is very sen- 
sitive to the r-min assumed. For example, if rmin ~ 0.3 
Mpc/h, the predicted virialized gas fraction would 
be fh ~ 50% and thus contradicts with our simula- 
tion result. One can also compare the prediction of 
p{> T) of halo model with our simulation. The gas 
temperature T of a halo mass m can be estimated 
by T = T»(m/miif^-\ where at z = 0, Ts ~ 5S1„ ~ 1.3 
keV (Pen.l998i^ and mg is the mean mass contained 
in a sphere with radius 8 Mpc/h. Thus the halo 
model predicts p(> l.SkeV) = fhija ^ mg) ~ 4%, and 
p(> O.lkcV) ~ fh{m ^ 0.02m8) ~ 50%. These predic- 
tions are roughly consistent with our simulation re- 
sult. Thus, a low fraction of virialized gas may be 
caused by simulation resolution. 

The effect of resolution to the p — T relation 
can be estimated as follow.^. The density profile of 
each halo can be approximated as S = Svii{r/rvii ~ 
100(r/rvir). The smoothing caused by limited resolu- 
tion convolves this density field with a window func- 
tion with radius Rf. For halos with rvir ^ Rf, the 
smoothed density 5^' will be 5^' ~ 100{r^ir/Rff oc T. 
Here, we have utilized that rvir oc 
For WMAP cosmology at 2 = 0, this implies that 
S ~ 100(lMpc/h/7?/)2(r/keV). This relation is sur- 
prisingly consistent with our simulation, if taking 
Rf — 0.2 Mpc/h. So the p — T relation may also be 
caused by artifacts of simulations and has to be used 
in caution. 

On the other hand, the simulation resolution is 
hard to explain all aspects of simulation results. The 
p-T relation shows a robust convergence in the self 
similar simulations (e.g. 20%, 50% p-T contour in 
Fig. 4) and only changes very slowly as the non- 
linear length scale changes by a factor of 4. This 
suggests that it is unlikely that the observed p-T re- 
lation is purely caused by simulation artifacts, oth- 
erwise the p-T relation would change significantly 
with resolution. 

These issues deserves further investigations. But 
more careful comparison requires more realistic halo 
mass function than the Press-Schechter formalism 
and a series of simulations to have fair sample of ^ 
keV halos and thus beyonds the scope of this paper. 

But despite these significant resolution effects on both 
density and temperature fields, the simulated overall p{> T) 
shows a good convergence except at very low or high tem- 
perature ranges (Fig. |5]&|SJ. Since the energy conservation 
guarantees the total amount of kinetic and thermal energy 
to be well simulated, the good agreement of p(> T) implies 
that the conversion efficiency from kinetic energy to thermal 
energy is well simulated too. Since Tg = J Tdp, the effect 



^ The discussion presented here is based on the comments of the 
anonymous referee. 
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Figure 6. The overall temperature distribution p(> T) for 
n = —I self similar simulation. For the temperature field, the 
simulation effect is minor except for very low or high temperature 
regions. Since the gas density weighted temperature fg = J Tdp, 
Tg is well simulated. 

of simulation limitations to Tg should be minor, as we will 
quantify in the next section. 

The simulated Tg (z) would be less resolution dependent 
and robust. We will develop our model for Tg{z) in the next 
section, test it against simulations, quantify simulation ar- 
tifacts and provide a precision model of Tg(z). As we will 
show in the simulation effect on Tg is only non-negligible 
at z < 0.5 and this effect is no bigger than ~ 20%, even at 
z = 0. 



4 THE CONTINUUM FIELD MODEL 

In a gravitational heating scenario, the gas temperature is 
determined by the gravitational potential $. The pressure 
depends on the thermalized fraction of the total kinetic en- 
ergy. The translational kinetic energy is thermalized from 
the energy released when particles shell cross. A model of 
the thermalized energy is thus given by the difference in en- 
ergy between two particles separated by a non-linear scale in 
Lagrangian space, which is the distance at which they can be 
expected to have shell crossed. The exact procedure amounts 
to solving the non-linear evolution equations directly. But we 
can treat the effect statistically in a linear fashion. In the 
initial linear evolution, the gravitational potential remains 
constant. After virialization, the gravitational energy at a 
fixed location remains almost constant. In an Eulerian de- 
scription, we can describe the energy of particles at a final 
virialized location as the energy released as a particle travels 
from its initial position to the final virialized location. We 
can then relate the gas temperature to $ through the viral 
theorem: 




1 + z 

Figure 7. The density weighted gas temperature Tg. Data points 
are the results of our self similar simulations (512^), WMAP 
simulation (1024^, 100 Mpc/h) and solid lines are our model 
predictions. Our model predicts a self similar scaling relation 
fg(z) oc (1 -I- z)("~l)/("+3) for self similar simulations. Our pre- 
dictions agree with simulations very well at most redshifts. The 
breaking of the self similar relation at low redshifts suggests the 
limitation of simulations. We will show that the discrepancy of 
the WMAP simulation at low redshifts is caused by simulation 
resolution in fig.|S]while the discrepancy of the self similar simual- 
tions is caused by limited box size. 

fcsT.(x) = [*(x) - $(x)] . (6) 

Since the initial position is not exactly known, we take a 
spherical average over the non-linear scale to average over 
all possible initial locations and obtain the mean initial po- 
tential 

$(X) =y"$(x')We(x-x')d'^2:'- (7) 

For a detailed explanation, refer to IZhang fc PenI i200ll) . 

Then, the averaged gas density weighted temperature 

Aig(fc)/e(fc)— . (8) 

Here, f^{k) = [1 - We{k)]/k^ [k is in unit of Mpc/h). A^g 
is the dark matter-gas cross correlation power spectrum. 

In our model, We{k) is a free function, but its asymp- 
totic behavior toward = is fixed by the requirement that 
Tg follows the density field at large scales, or equivalently, 
the Tg bias with respect to the underlying density field is 
a constant at large scales. Its behavior at small scales is 
hard to determine from first principles. But since at scales 
smaller than smoothing scale, We{k) — » and fe{k) —* fc~^, 
the exact behavior of We{k) at large k is not very important. 
Based on these considerations, a natural choice of We{k) is 
a Gaussian function We{k) = exp(— fc^r^). For this func- 
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Figure 8. The resolution effect of simulations to Tg. The data 
points with {2a) error bars are the WMAP simulation result while 
the triangle data points are our model prediction using the simu- 
lated . For clarity, triangle data points are shifted horizontally 
arbitrarily. The solid line is our model prediction assuming gas 
perfectly follows dark matter and the dash line is the prediction 
assuming a small scale cutoff in the gas density power spectrum. 
The excellent agreement between two sets of data points supports 
the validity of our model and implies the simulation limitations 
to be the cause of the apparent discrepancy found in fig. |3 The 
excellent reproduction of the simulation results at low redshifts 
by the dash line implies that the simulation resolution is the cause 
of the discrepancy. 



tion, when fc 0, fe{k) — » r^, so the temperature bias with 
respect to the density field is a constant. Since gas gains 
thermal energy by shell crossing, which happens at the non- 
linear scales, we expect to be roughly equal to the density 
correlation length. So, the evolution of re{z) follows that of 
the density correlation length. Essentially, the only free pa- 
rameter in our model is re{z = 0). We will choose = tnl, 
where rjvL is the non-linear scale set by ^l(?'nl) = 1 with 
as the linear correlation function. This choice of re has to 
be tested against simulations and this is the only parameter 
that requires calibration against simulations. 

For self-similar simulations, A5g(fc) should scales as 
/(fc/^NL) (We define A;nl as /cnl = I/^nl), then Eq.|Hl nat- 
urally predicts 

fg{z) = foiz = 0)(1 + (9) 

Without feedback or cooling, the gas should follow the dark 
matter distribution matter to very high overdensity and thus 
we expect A'i^ik) = Al^(k) = A^(A;). 

W e calculate A^^{k) using the code of ISmith et al. I 
l|2003l) . We compared the predictions from our model with 
simulation results and found a good agreement (Fig.|7J. The 
scaling relation Eq. (|5J with the right amplitude is observed 



Figure 9. The gas-dark matter relation. The dash lines are the 
gas-dark matter cross correlation coefficients while the solid lines 
are the gas biases versus dark matter at different redshifts. At 
small nonlinear scales, gas ceases to follow dark matter in simula- 
tions. But the redshift dependence of such behavior suggests that 
it is unphysical and possibly caused by the simulation resolution 
effect. Lower redshift simulation has better mass resolution (left 
panel), which corresponds to better spacial resolution in unit of 
nonlinear scale. But spacial resolution in unit of absolute physical 
scale is poorer at lower z (right panel). 



at high redshifts and further confirms the validity of our 
model. 

At low redshifts, the scaling relation breaks down for 
the self similar simulations. This is caused by the finite sim- 
ulation box size. Its effect to Tg corresponds to a lower k 
cutoff fccut = 2n/L in the integral of Eq. |H] where L is the 
box size. For our self similar simulations, the correlation 
length at z = is half the box size and /cnl = 2/L. One has 
fccut ~ vrfcNL- So, it is the limited box size that causes Tg in 
self similar simulations to lose power at low redshifts. 

The deviation between the predicted and simulated Tg 
is also observed in the WMAP simulation. Such discrepancy 
increases toward low redshifts and exceeds 2a level at z = 
(Fig. |HJ , so it is hard to be explained by sample variance. 
If this discrepancy is caused by simulation limitations, Tg 
calculated using the simulated WMAP A^g(fc) should agree 
with the simulated Tg. Indeed, the agreement is better than 
5% at low redshifts (Fig. (HJ. Thus we show that this dis- 
crepancy can be naturally explained by the simulation limi- 
tations and thus our model works well to better than several 
percent at low redshifts. 

We further probe which simulation limitation causes 
this discrepancy. For WMAP simulations, even at z = 0, 
the nonlinear scale is still much smaller than the box size, 
so the box size effect is negligible. Resolution effect causes 
A^g to lose power at small scales and causes the simulated 
Tg to lose power. Since the resolution of the hydro part of 
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a hydro simulation is generally worse than that of the N- 
body part, gas ceases to follow dark matter below certain 
scale. Such deviation is a suitable measure of simulation 
resolution and can be quantified. We define the gas bias 
bg{k) = Ag (k) / A'^^_^{k) and the gas-dark matter cross 
correlation coefficient r = A|g(fc)/iyA|^^^"(fc)Af(fcy. We ex- 
pects bg{k) < 1 at very nonlinear scales. This behavior is 
observed in our simulations (Fig. O. Poorer resolution of 
gas with respect to dark matter means gas is smoother than 
dark matter, so phenomenologically, one can treat the gas 
density field as a smoothing of the underlying dark matter 
density field: 

Sg{x)^ J 5in.{^')Wg{X-^')d^x' . (10) 

In this model, gas perfectly correlates with dark matter and 
r = 1. In our simulation, we find that at A^^^ < 200, this 
is the case (Fig. [SJ. One can model Wg{k) — exp[— fc^/fcg]. 
An ideal simulation should have a kg such that A^{kg) ^ 1. 
In simulations, kg{z) should increase with z since for higher 
z, the nonlinear scale is smaller (e.g. fig. The simulated 
kg can be modeled by kg{z) — 5(1 -I- z)^ h/Mpc, which is 
roughly consistent with the simulated gas power spectrum, 
and reproduces the simulation results (Fig. (HJ. This agree- 
ment implies that for WMAP cosmology, the simulation res- 
olution causes the simulated Tg to lose power at low redshift. 

This conclusion seems counter-intuitive. Since gas tem- 
perature arises from thermalization at nonlinear scales, 
which are better resolved in lower redshift simulations (left 
panels of fig.|5]and bottom panel of fig. llUL we may expect 
less severe resolution problem for simulated Tg at lower red- 
shifts. For self similar simulations, this is true since Ainl is 
the only relevant scale. But ACDM cosmology breaks self 
similar condition and makes the nonlinear scale ^nl not the 
only relevant parameter to determine Tg . We show how this 
running power index causes more severe resolution problem 
for simulated Tg. 

The relative contribution from different scale k to Tg re- 
lies on the slope of the power spectrum. The larger the power 
index at fc > /cnl is, the larger the relative contribution to 
Tg is and the higher the requirement of spacial resolution in 
unit of /cnl to simulate Tg is. For WMAP cosmology, /cnl 
keeps decreasing toward low z. The effective power index 
at fc = /cnl keeps increasing and the relative contribution 
from k > fcNL keeps increasing. Such behavior requires a 
stronger spacial resolution in unit of ^nl. In order for the 
simulated Tg not to lose power, at z — 1, simulation must 
be able to resolve k < 2feNL, but ai z = 0, the requirement 
is fc < IO/cnl (Fig. 1101 . As we see from the bottom panel of 
Fig. 1101 WMAP simulation meets this requirement ai z = 1 
but fails at 2; = 0. So, we conclude that, for CDM simula- 
tion with reasonable large simulation box (> 50 Mpc/h), the 
simulation resolution is the dominant simulation limitation 
to simulate Tg. 

In the above discussion, we have assumed that 
in the relevant k range, Ag — A^^ = A^g. Though 
at highly nonlinear regime, the simulated Ag does 
lose power with respect to A'^^ (Fig. this is 
likely caused by resolution effect since when increas- 
ing resolution, = Ag/A^^ keeps increasing toward 
unity (Fig. 9). Since we do not have higher resolution 
simulations to test if b will reach unity in these non- 
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Figure 10. The requirement of Tg on simulation resolution. The 
top panel plots the integrand of Eq.|S] namely, the contribution of 
differen t scales to Tg . The n onlinear power spectra are calculated 
by the ISmith et aU <2003l) code. At z > 1, only k ~ /cnl is 
required to resolve. But at 2: = 0, A; > IO/cnl is required to resolve. 
The bottom panel is the gas bias in our simulation. Its deviation 
from unity is a measure of the simulation resolution. We find that, 
at 2 > 1, WMAP simulation meets the resolution requirement. 
But at 2 = 0, the resolution requirement [k ~ IO/cnl) is beyond 
simulation ability. Since simulated power spectrum loses power 
at 10% level at fc ~ IO^nl, we expect the simulated Tg to lose 
power at 10% level, as predicted in Fig. ISl 



linear regimes, this question still remains open. But 
one can estimate it using the halo model. One ex- 
pects that the significant deviation from 6=1 could 
only happen where gas pressure is comparable to 
gravity, namely, in each halos. Observationally, intr- 
acluster gas develops a constant density core while 
dark matter may develop a density cusp toward the 
center. It is not clear whether gas pressure alone 
could generate such gas core or other mechanism 
such as feedback is required. But even the devia- 
tion in the gas and dark matter distribution presents 
in adiabatic simulations (with infinite resolution), it 
should only happen at very small scale r < rc ~ 0.1 
Mpc/h, or fe > 10/i/Mpc. Thus we expect that 6 = 1 is 
a good approximation at scales k < 10/t/Mpc. Since 
the dominant contribution to Tg for WMAP cosmol- 
ogy comes from k <C 10/i/Mpc, we do not expect the 
possible real deviation of b to be responsible for the 
discrepancy between predicted and simulated Tg. We 
thus neglect this possible effect. 

In summary, our choice of re, the only free parameter 
in our model passes the tests of all simulations. Thus our 
model has no free parameter, is free of simulation artifacts 
and is able to predict the real Tg to several percent accuracy. 
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Figure 11. The predicted Tg and y as a function of erg and 
Sim- For the left panels, from top down, lines correspond to 
f2m = 0.5,0.4,0.3,0.2. For the right panels, from top down, lines 
correspond to erg = 1-1, 1-0, 0.9, 0.8, 0.7. 



It is interesting to compare these results with 
the halo model prediction. For the self similar 
case, the same scaling relation (Eq. ^ is predicted. 
The scaling relation does not depend on the spe- 
cific form of the halo mass function, either Press- 
Schechter formalism or the Jenkins fitting formula. 
It only relies on (1) the virial theorem Tg oc m^/rvir 
and (2) the general form of halo mass function 
dn/dm oc m^'^ f(iy')\dlnu/dlnm\. The detailed study 
(Refreeier &: Tevssier 2002) shows that the halo 
model prediction is consistent with simulations. Sur- 
prisingly, the halo model also tends to predict a 
higher Tg at low z than simulations and strongly 
implies the same resolution eff'ect. The continuum 
model and the halo model rely on different as- 
sumptions and require different input. Though both 
models are based on virial theorem, the continuum 
model treats gas as continuum field and the halo 
model treats gas as distributed in discrete halos. 
The continuum model requires the nonlinear den- 
sity power spectrum, while the halo model requires 
the halo mass function and the halo m — T relation. 
The consistency in the predictions of two distinc- 
tive models gives us confidence that the predicted 
Tg is reliable and simulations artifacts are suitably 
handled. 

For the WMAP cosmology, we predict fg{z = 0) = 0.32 
keV. Tg is sensitive to erg and Q,rn- For self similar cosmology, 
Tg OC cr-("-i'''("+^). For ACDM, the actual dependences of 
Tg on (jg and Q,m (fig. lllH is complicated due to the running 
index of the density power spectrum. In our interested erg 
and Sim range, the effective power index licit is —2 < ricff < 



— 1. So Tg oc (Tg with a ~ — (rzeff — 1)/ (jieB+S) ~ 1-3 and 
/3 ~ 1. A smaller erg results in a larger Anl and thus a smaller 
effective power index n^f!. So a is larger. The deviation of /3 
from unity comes from the dependence of Wcff on fim since 
CDM transfer function depends on the combination q ~ 
k/Q.m- Around the WMAP cosmology Q,m ~ 0.268 and as = 
0.84, Tg can be fitted as 

Tg = O.S2{as/0.84:f-°^~°-^^^'"{n^/0.268y-'^^-°-'^'''' keV.(ll) 



5 THE SZ MEAN Y PARAMETER 

The SZ mean y parameter is calculated using Eq. |21 and 
the result is shown in Fig. 1111 y is generally ~ 10~^. For 
the WMAP cosmology, y — 2.6 x 10~® and the mean tem- 
perature decrement at Rayleigh- Jeans regime is 14/iK. The 
dominant contribution comes from z ~ 1 fFig. I12II . Since at 
2 > 1/2, dx/dz oc l/\/flrn, oue may expect y oc fg/y/Tl^ oc 
But since in a higher matter density universe, the 
density field evolves faster and thus Tg drops faster with in- 
creasing z, the y dependence on Qm is stronger than Q^'^. 
Indeed, we find y oc Sl^^. Around the WMAP cosmology 
f^m ~ 0.268 and erg = 0.84, y can be fitted as 

y = 2.6 X 10"'^(erg/0.84)''-'"^"" (nm/0.268)'-''**-'^-'^"** . (12) 

Though our model is able to predict the Compton y 
in an adiabatically evolving universe to several percent ac- 
curacy, it does not include any non-gravitational thermal 
processes, which introduce non-negligible effect to y. Photo- 
ionization contributes j/photonion ~ rlO^K/meC^ ~ 3 x 10~^, 
or ~ 10% of the adiabatic IGM y. Though feedback, pre- 
heating, radiative cooli ng may decrease the SZ power spec- 
trum by a factor of 2 Jda Silva et al.ll200lt iLin et aljl2002l : 
IWhite. Hernauist fc SpingeM2 002': 'Zhang fc Wu"2003^. they 
only affects y at 10% level fe.g. White. Hernauist fc Sp i ngell 
(2002)). This is straightforward to understand. Once hydro- 
static equilibrium is reached, the gas pressure is always de- 
termined by the gravitational potential, which is mainly set 
by dark matter distribution and is only weakly affected by 
these thermal processes. Feedback and preheating do not 
change the total amount of gas. While radiative cooling 
turns some gas into bound objects, such mass loss is minor 
since bounds objects in galaxies only accounts for < 10% 
baryons. Thus, the change of the total gas thermal energy 
due to these processes is minor. These processes change y 
mainly during the stages of expansion in the case of feed- 
back and preheating and infall in the case of radiative cool- 
ing. Such stages are either in semi hydrostatic equilibrium 
(feedback and radiative cooling) or last relatively short time 
(mild preheating), thus their effects to y is not sign ificant. 
Clust er magnetic field also only has ~ 10% effect to y llZhand 
I2004D . 

WMAP measured a high Thomson optical depth to 
the last scatter surface and implies an early reionization 
epoch caused by first stars. At high z, CMB density is 
high and is able to convert a considerable fraction of first 
star supernova explosion thermal energy through the ef- 
ficient Compton cooling. Such first star contribution to 
r/ is ~ few 10~^ and comparable to low redshift IGM y 
iOh. Coorav fc Kamionkowskill2003f) . 

These processes have distinctive signatures to y and the 
SZ power spectrum, respectively. For first stars, since at high 
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Figure 12. The integrand of y, as defined by dy/dlnz for the 
WMAP cosmology. The dominant contribution to y comes from 
2 ~ 1. 

z ~ 10, density fluctuations are small, the relative contribu- 
tion of first stars to the SZ fi uctuation is much smaller than 
that o f to y. As estimated bv lOh. Coorav &: Kamionkowskil 
i2003h . even if first star y is larger than that of low redshift 
IGM, its contribution to the SZ power spectrum can still be 
an order of magnitude smal ler than that of low redshift IGM 
(Fig. 1, lOh, Coorav fc Kamionkowskil feOOgh l. The case of 
photon-ionization is similar, but that of feedback, preheat- 
ing, radiative cooling and magnetic field is opposite. Thus, 
combining the y and the power spectrum measurement helps 
to separate these contributions more unambiguously. For ex- 
ample, if y can be measured to ~ 10% accuracy, the first 
star contribution can be constrained with a statistical un- 
certainty ~ 0.4 X 10~® and systematic underestimation of 
~ 0.5) X 10~® caused by feedback, etc. 

Unfortunately, the direct precision measurement of ab- 
solute y is difficult. Currently, the best measurement, y < 
1.5 X 10~^ is given by the COBE/FIRAS measurement 
iFixsen. et alj[l996l) . This result has already been able to 
rule out combinations of high as > 1.1 and high Urn 0.5. 
This constrain is quite weak. But considering the con- 
tributions from non-gravitational processes could make it 
stronger. Nonetheless, y may be measured to a higher accu- 
racy in the future and/or inferred from new statistics and 
helps to independently constrain as and $lm . 



tics. The two ways of the SZ modeling, analytical models 
and hydro simulations both have their own limitations. It is 
essential to quantify simulations limitations and test analyt- 
ical models against corrected simulations. The convergence 
of Tg stands as the lowest requirement for simulations to 
reliably predict the SZ effect. 

We ran n = —1,-2 self similar 512^ hydro simulations 
to quantify simulation limitations utilizing their self similar 
scaling relation. We also ran a high resolution 1024^ cell, 
100 Mpc/h box size hydro simulation adopting WMAP cos- 
mology. We find that the simulated p(> T), the fraction of 
mass with temperature bigger than T, shows a good con- 
vergence for all our simulations, except at both tails. This 
convergence suggests that Tg is well simulated. Our contin- 
uum field model is then tested against these simulations. It 
passed all tests and we believe that its prediction for y is 
accurate to several percent. 

Various simulation limitations such as limited box size 
and limited resolution can affect the simulated Tg. Gener- 
ally, for a ACDM simulation, the nonlinear scale is much 
smaller than the box size, thus the box size effect is neg- 
ligible and the resolution effect is the dominant cause of 
simulation artifacts in Tg. Due to curvature in the CDM 
power spectrum, the temperature is more difficult to model 
accurately in simulations at fixed comoving resolution. We 
found that, at z = 0, due to the simulation resolution, gas 
power spectrum loses power at small scales with respect to 
dark matter power spectrum. This behavior causes the sim- 
ulated gas density weighted temperature Tg to be ~ 20% 
underestimated. But this resolution effect becomes negligi- 
ble quickly toward higher redshift. At z > 0.5, simulated 
Tg is quite accurate. Since the dominant contribution to y 
comes from 2 ~ 1, our simulation prediction of y is reliable 
to ~ 10% level. Furthermore, our analytical model is able to 
correct this simulation artifacts and predicts y with several 
percent accuracy. For a fiat, low matter density ACDM uni- 
verse, y = 2.6 X 10-'5(cr8/0.84)'' i-2"'"(n™/0.268)i-2*-°-2"8. 
The current upper limit ofy < 1.5 x 10~^ measured by FI- 
RAS has already ruled out combinations of high erg > 1.1 
and high Qm ^ 0.5. 

Our simulations confirms previously found T oc p re- 
lation in a large region of p-T plane. Simulation resolution 
may be partly responsible for this relation. But the results 
of self similar simulations show robust convergence of this 
relation and can not be explained by simulation resolution 
artifacts. This issue deserves a further investigation. We also 
found that, the simulated p{> T,> 6), does not converge. At 
high density regions, it is caused by density resolution lim- 
itation while at low density regions, it is caused by failure 
of capturing shocks. Though this numerical limitation has 
only minor effect on Tg, it may affect the gas pressure power 
spectrum a lot. This issue will be addressed in a companion 
paper. 



6 SUMMARY 

The mean gas density weighted temperature Tg and the 
mean SZ Compton y parameter are the lowest order SZ 
statistics. Their precision modeling stands as the first step 
toward the precision understanding of the IGM SZ effect and 
may provide useful clues for modeling of higher order statis- 
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